import os
#from photometry import *
import numpy
from optparse import OptionParser
import os.path
import glob
from types import NoneType
import matplotlib.pyplot as plt
from chimera.core.chimeraobject import ChimeraObject
import datetime as dt

__author__="gilfus"
__date__ ="$Jul 8, 2009 8:54:57 PM$"
VERBOSE = True

class Extinction(ChimeraObject):
    __config__ = {'database': '/Database/0',
			 'filters': 'U,B,V,R,I',
			 'referenceBand': 'V'}
    magnitudes = []
    airmasses = []
    colorIndices = []
    m0 = []
    a = 0
    k = 0
    kprime = 0
    def __init__(self):
        ChimeraObject.__init__(self)

    def __start__(self):
        self.database = self.getManager().getProxy(self['database'])

    #For a given night and filter, calculate a and k for a curve of best fit of
    #    the form m_inst = m0 + a + k*AM,
    #    where m0 for each star is known and is retrieved from the database.
    def getExtCatalog1stOrder(self, filterIn, midnight, catalog, plotFile=None,
					amLbound=None, amUbound=None,
					magLbound=None, magUbound=None,
					compareMags=False):
		#Clear self.magnitudes, self.airmasses
		self.magnitudes = []
		self.airmasses = []
		
		#Compute night start time and night end time for night during which 
		#'time' falls
		nightStartTime = midnight - dt.timedelta(hours=12)
		nightEndTime = midnight + dt.timedelta(hours=12)
		
		#Make a list, expStars, of all catalog star observations made on the 
		#given night through the given filter
		catalogStars = self.database.getStars(catalog=catalog, 
								   hasExposures=True, 
								   magFilter=filterIn)
		if len(catalogStars) == 0:
			return ['Error: No catalog stars found tonight.  '+
				   'Cannot compute extinction.']
		expStars = []
		magLB = None
		magUB = None
		for star in catalogStars:
			if magLbound!=None: magLB = magLbound+star.magnitudes[filterIn]
			if magUbound!=None: magUB = magUbound+star.magnitudes[filterIn]
			expStars.extend(self.database.getStarExposures(star=star,
					    filter=filterIn, startObsDate=nightStartTime, 
					    endObsDate=nightEndTime, minAM=amLbound, maxAM=amUbound,
					    minMag=magLB,
					    maxMag=magUB))
            
		#For each star in expStars, append its airmass to self.airmasses and
		#its m_inst - m0 to self.magnitudes
		obsMag = list()
		zeroMag = list()
		for n in range(len(expStars)):
			self.airmasses.append(float(expStars[n].airmass))
			obsMag.append(float(expStars[n].magnitude))
			zeroMag.append(float(self.database.getStarMagnitude(expStars[n], filterIn)))
			self.magnitudes.append(obsMag[-1] - zeroMag[-1])
		
		if compareMags:
			self.compareMagnitudes(expStars=expStars)
		
		#Use numpy.polyfit to calculate a and k
		self.k = numpy.polyfit(self.airmasses, self.magnitudes, 1)[0]
		self.a = numpy.polyfit(self.airmasses, self.magnitudes, 1)[1]
		
		#Calculate correlation coefficient
		ssxx = 0
		ssyy = 0
		ssxy = 0
		xmean = 0
		ymean = 0
		for airmass in self.airmasses:
			xmean += airmass
			ymean += airmass
		xmean /= len(self.airmasses)
		ymean /= len(self.magnitudes)
		for n in range(len(self.airmasses)):
			ssxx += (self.airmasses[n] - xmean)**2
			ssyy += (self.airmasses[n] - ymean)**2
			ssxy += (self.airmasses[n] - xmean) * (self.magnitudes[n] - ymean)
		r2 = (ssxy**2)/(ssxx*ssyy)
		
		if isinstance(plotFile, str):
			plt.clf()			
			plt.suptitle('1st Order Extinction: m-m0 = '+str(self.a)+' + '+str(self.k)+'AM')
			plt.xlabel('airmass')
			plt.ylabel('m-m0')
			plt.plot(self.airmasses,self.magnitudes,'b.')
			plt.plot([min(self.airmasses)-.02, max(self.airmasses)+.02], 
				   [self.a+self.k*(min(self.airmasses)-.02), 
				    self.a+self.k*(max(self.airmasses)+.02)], 'r-')
			plt.savefig(plotFile+'1stO.png')
		
		return {'k':self.k, 'a':self.a, 'correlation':r2}
	
    #For a given night and filter, compute coefficients a k and k' for a surface 
    #of best fit of the form m_inst = m0 + a + k*AM + k'*AM*color
    #color is the difference between filterIn and the reference band (i.e. V for 
    #   the Johnson system)
    #   The redder band is subtracted from the bluer band (i.e. color = U-V or 
    #		color = V-R)
    #   If filterIn is the reference band, the next-bluer band is used (if 
    #		filterIn = V, color = B-V)
    #Saves two plots to plotFile1stO.png and plotFile2ndO.png
    def getExtCatalog2ndOrder(self, filterIn, midnight, catalog, plotFile=None, 
		amLbound=None, amUbound=None, magLbound=None, magUbound=None,
		colorLbound=None, colorUbound=None, compareMags=False):
		#clear self.magnitudes, self.airmasses, self.colorIndices
		self.magnitudes = []
		self.airmasses = []
		self.colorIndices = []
		
		#Compute night start time and night end time for night during which
		#'time' falls
		nightStartTime = midnight - dt.timedelta(hours=12)
		nightEndTime = midnight + dt.timedelta(hours=12)
		
		#Get a list of the filters used tonight and find the index of the 
		#reference band in this list
		referenceBand = self['referenceBand']
		bands = self['filters'].split(',')
		filterList = list()
		hasRefBand = False
		for band in bands:
			filterList.append(band)
			if band == referenceBand: hasRefBand = True
		if not hasRefBand: return ['Filter system does not contain '+
							'reference band \''+referenceBand+'\'.']
		if len(filterList) <= 2:
			return ['Only one filter in system.  To calculate extinction '+
					'system must contain at least the reference band '+
					'and the input filter.']
		referenceIndex = -1
		for n in range(len(filterList)):
			if filterList[n] == referenceBand:
				referenceIndex = n
		if referenceIndex == -1:
			return ['Error: reference band, \''+referenceBand+'\', not '+ 
				  'found in list of tonight\'s filters.  Cannot '+
				  'compute extinction.  Tonight\'s filters: '+filterList]
		if referenceIndex == 0:
			return ['Error: reference band, \''+referenceBand+'\', bluer'+
				  ' than all other filters used tonight.  Cannot '+
				  'compute extinction.  Tonight\'s filters: '+filters]
		if filterList.index(filterIn) < referenceIndex:
			ciString = filterIn+'-'+referenceBand
		elif filterList.index(filterIn) > referenceIndex:
			ciString = referenceBand+'-'+filterIn
		else:
			ciString = filterList[referenceIndex-1]+'-'+referenceBand
		
		#Make a list, expStars[], of all catalog star observations made on
 		#the given night through the given filter
		catalogStars = self.database.getStars(catalog=catalog, 
								   hasExposures=True,
								   magFilter=filterIn)
		#Remove stars with no exposures in filterIn tonight or without
		#necessary m0s
		n = len(catalogStars) - 1
		while n >= 0:
			if len(self.database.getStarExposures(
								star=catalogStars[n],
								filter=filterIn,
								startObsDate=nightStartTime,
								endObsDate=nightEndTime))==0:
				catalogStars.pop(n)
			elif filterIn != referenceBand:
				if self.database.getStarMagnitude(catalogStars[n],
									     referenceBand)==None:
					catalogStars.pop(n)
			else:
				if self.database.getStarMagnitude(catalogStars[n],
								filterList[referenceIndex-1])==None:
					catalogStars.pop(n)
			n -= 1
		if len(catalogStars) == 0:
			return ['Error: No stars in catalogue \''+catalog+'\' found '+
				  'tonight.  Cannot compute extinction.']
		expStars = []
		starExposures = [] #Number of ExposureStar s for each star in 
					 #catalogStars; only used if plotFile!=None
		
		#Find out how many ExposureStars will not be used for each Star
		magLB = None
		magUB = None
		for star in catalogStars:
			esLen = len(expStars)
			if magLbound!=None: magLB = magLbound+star.magnitudes[filterIn]
			if magUbound!=None: magUB = magUbound+star.magnitudes[filterIn]
			expStars.extend(self.database.getStarExposures(star=star,
					filter=filterIn,
					startObsDate=nightStartTime,
					endObsDate=nightEndTime, minAM=amLbound,
					maxAM=amUbound,
					minMag=magLB,
					maxMag=magUB))
			if isinstance(plotFile, str):
				starExposures.append(
					len(self.database.getStarExposures(star=star,
					filter=filterIn, startObsDate=nightStartTime,
					endObsDate=nightEndTime, minAM=amLbound,
					maxAM=amUbound, minMag=magLB,
					maxMag=magUB)))
		if len(expStars)==0:
			return ['Error: No stars in catalogue \''+catalog+
				  '\' have images in band \''+filterIn+'\'.']
		multipleCI = 0
		
		#For each exposure star in expStars, add its airmass to
		#self.airmasses, its color index to self.colorIndices, its magnitude 
		#to self.magnitudes
		zeroMag = list()
		obsMag = list()
		for n in range(len(expStars)):
			self.airmasses.append(float(expStars[n].airmass))
			if filterList.index(filterIn) < referenceIndex:
				self.colorIndices.append(
					float(self.database.getStarMagnitude(
					expStars[n].star, filterIn) -
					self.database.getStarMagnitude(expStars[n].star,
										 referenceBand)))
			elif filterList.index(filterIn) > referenceIndex:
				self.colorIndices.append(float(self.database.getStarMagnitude(
					expStars[n].star, referenceBand) - 
					self.database.getStarMagnitude(expStars[n].star, filterIn)))
			else:
				self.colorIndices.append(float(self.database.getStarMagnitude(
					expStars[n].star, filterList[filterList.index(referenceBand)-1]) - 
					self.database.getStarMagnitude(expStars[n].star, referenceBand)))
			#Check that we have observations of at least two different color indices
			if n > 0 and self.colorIndices[n] - self.colorIndices[n-1] != 0:
				multipleCI = 1
			obsMag.append(float(expStars[n].magnitude))
			zeroMag.append(float(self.database.getStarMagnitude(expStars[n].star, filterIn)))
			self.magnitudes.append(obsMag[-1] - zeroMag[-1])
		if multipleCI == 0:
			return ['All stars in catalog \''+catalog+'\' observed tonight have color index'+
				  '('+ciString+') = '+str(self.colorIndices[0])+'.  '+
				  'Must have at least two different color indices to compute extinction.']
		
		if compareMags:
			self.compareMagnitudes(expStars=expStars)
		
		#Compute a, k, and k' using method of least squares
		#1.  Write RSS = Sum[(z - a - k*x - k'*x*y)**2], x=airmass, y=colorindex, z=m_inst-m0
		#2.  Set dRSS/da = 0, dRSS/dk = 0, dRSS/dk' = 0
		#    Now we have 3 linear equations and 3 unknowns:
		#        A1a + B1k + C1k' + D1 = 0
		#        A2a + B2k + C2k' + D2 = 0
		#        A3a + B3k + C3k' + D3 = 0
		#    Where all coefficients are in terms of observed x, y, and z
		#3.  Compute the sums needed to calculate coefficients
		N = len(self.magnitudes)
		xsum = 0.0
		x2sum = 0.0
		zsum = 0.0
		xysum = 0.0
		xzsum = 0.0
		x2ysum = 0.0
		xyzsum = 0.0
		x2y2sum = 0.0
		for n in range(N):
			xsum += self.airmasses[n]
			x2sum += self.airmasses[n]**2
			zsum += self.magnitudes[n]
			xysum += self.airmasses[n]*self.colorIndices[n]
			xzsum += self.airmasses[n]*self.magnitudes[n]
			x2ysum += (self.airmasses[n]**2)*self.colorIndices[n]
			xyzsum += self.airmasses[n]*self.colorIndices[n]*self.magnitudes[n]
			x2y2sum += (self.airmasses[n]**2)*(self.colorIndices[n]**2)
		#4.  Put coefficients into lists:
		#        row1 = [A1, B1, C1, D1]
		#        row1 = [A2, B2, C2, D2]
		#        row1 = [A3, B3, C3, D3]
		row1 = [-1.0, -1.0*xsum/N, -1.0*xysum/N, zsum/N]
		row2 = [-1.0*xsum/x2sum, -1.0, -1.0*x2ysum/x2sum, xzsum/x2sum]
		row3 = [-1.0*xysum/x2y2sum, -1.0*x2ysum/x2y2sum, -1.0, xyzsum/x2y2sum]
			#We have already checked that N!=0 and x2y2sum!=0
		#5.  Check that there is a unique non-trivial solution.  There are two cases:
		#5a. If the system of equations is homogeneous (D1=D2=D3), 
		#	then either there is an infinite solution set or the only solution is a=k=k'=0.
		#	Neither of these is acceptable.
		if row1[3]==row2[3] and row3[3]==row[2]:
			return ['Homogeneous system; no unique non-trivial solution.']
		#5b. If the system is inhomoeneous, there is only a solution if the determinant of 
		#	the matrix of the coefficients is not zero, i.e.
		#        /A1 B1 C1\
		#    det( A2 B2 C2 ) != 0
		#        \A3 B3 C3/
		if (row1[0]*(row2[1]*row3[2]-row2[2]*row3[1]) +
		   row1[1]*(row2[2]*row3[0]-row2[3]*row3[0]) +
		   row1[2]*(row2[0]*row3[1]-row2[1]*row3[0]) == 0):
			return ['No unique solution to system of equations.  Cannot compute extinction.']
		f1 = row1[0]
		for n in range(len(row1)):
			row1[n] /= f1
		f1 = row2[0]/row1[0]
		f2 = row3[0]/row1[0]
		for n in range(4):
			row2[n] -= row1[n]*f1
			row3[n] -= row1[n]*f2
		f1 = row2[1]
		for n in range(4):
			row2[n] /= f1
		f1 = row3[1]/row2[1]
		f2 = row1[1]/row2[1]
		for n in range(4):
			row3[n] -= row2[n]*f1
			row1[n] -= row2[n]*f2
		f1 = row3[2]
		for n in range(4):
			row3[n] /= f1
		f1 = row1[2]/row3[2]
		f2 = row2[2]/row3[2]
		for n in range(4):
			row1[n] -= row3[n]*f1
			row2[n] -= row3[n]*f2
		self.a = -1*row1[3]
		self.k = -1*row2[3]
		self.kprime = -1*row3[3]
		
		#Calculate correlation coefficient
		ssxx = 0
		ssyy = 0
		ssxy = 0
		xmean = 0
		ymean = 0
		for n in range(len(self.airmasses)):
			xmean += self.airmasses[n]
			ymean += self.magnitudes[n]
		xmean /= len(self.airmasses)
		ymean /= len(self.magnitudes)
		for n in range(len(self.airmasses)):
			ssxx += (self.airmasses[n] - xmean)**2
			ssyy += (self.airmasses[n] - ymean)**2
			ssxy += (self.airmasses[n] - xmean) * (self.magnitudes[n] - ymean)
		r2 = (ssxy**2)/(ssxx*ssyy)
		
		#Make plot
		if isinstance(plotFile, str):
			fig = plt.figure()
			fig.set_figheight(2.0*len(catalogStars))
			fig.subplots_adjust(hspace=.6)
			sp = []
			minAM = min(self.airmasses)
			maxAM = max(self.airmasses)
			amRange = maxAM-minAM
			minMag = min(self.magnitudes)
			maxMag = max(self.magnitudes)
			magRange = maxMag-minMag
			minx = minAM-amRange*.1
			maxx = maxAM+amRange*.1
			miny = minMag-magRange*.1
			maxy = maxMag+magRange*.1
			numcatalogs = len(catalogStars)
			o1 = numpy.polyfit(self.airmasses, self.magnitudes, 1)
			amAll = []
			magAll = []
			for i in range(numcatalogs):
				m = 0
				largestCI = -100
				for n in range(len(catalogStars)):
					if self.colorIndices[m] > largestCI:
						largestCI = self.colorIndices[m]
						lgCIn = n
						lgCIm = m
					m+=starExposures[n]
				am = []
				mag = []
				n = lgCIm
				while n < lgCIm + starExposures[lgCIn]:
					am.append(self.airmasses[n])
					mag.append(self.magnitudes[n])
					amAll.append(self.airmasses[n])
					magAll.append(self.magnitudes[n])
					n+=1
				sp.append(fig.add_subplot(numcatalogs,1,i+1, title=ciString+' = '+
					str(largestCI)+': m_inst - m0 = '+str(self.a)+' + '+
					str(self.k+self.kprime*largestCI)+'AM', ylabel='m-m0'))
				l1 = sp[-1].plot(am, mag, 'bo')
				l2 = sp[-1].plot([minx,maxx],
						[self.a+self.k*minx+self.kprime*minx*largestCI,
						 self.a+self.k*maxx+self.kprime*maxx*largestCI], 'r-')
				l3 = sp[-1].plot([minx,maxx],[o1[1]+o1[0]*minx,o1[1]+o1[0]*maxx], 'g-')
				sp[-1].set_xlim(minx, maxx)
				sp[-1].set_ylim(miny, maxy)
				if i == numcatalogs - 1:
					sp[-1].set_xlabel('airmass')
			
				n = lgCIm + starExposures[lgCIn] - 1
				while n >= lgCIm:
					self.airmasses.pop(n)
					self.magnitudes.pop(n)
					self.colorIndices.pop(n)
					n -= 1
				catalogStars.pop(lgCIn)
				starExposures.pop(lgCIn)
			fig.legend((l1,l2,l3),('Observation', '2nd Order Prediction',
						     '1st Order Prediction'),'upper left')
			plt.savefig(plotFile+'2ndO.png')
			
			#Make a plot of AM v. m_inst-m0 of all stars
			fig = plt.figure()
			fig.suptitle('1st Order Extinction: m-m0 = '+str(o1[1])+' + '+str(o1[0])+'AM')
			sp2 = fig.add_subplot(1,1,1)
			sp2.plot(amAll, magAll, 'bo')
			sp2.plot([minx,maxx],[o1[1]+o1[0]*minx,o1[1]+o1[0]*maxx], 'g-')
			sp2.set_ylabel('m-m0')
			sp2.set_xlabel('airmass')
			plt.savefig(plotFile+'1stO.png')
		
		return {'a': self.a, 'k': self.k, 'kprime': self.kprime, 'correlation':r2}
    
    #For i stars whose established zero-airmass magnitudes are not known,
    #    do a least squares fit such that m_i = m0_i + k*airmass,
    #    where k is the same for all stars and m0_i is different for each star
    #Parameter starFile is the name of a file whose each line contains ra and dec of
    #    a star separated by a space
    #Saves two plots to plotFileall.png and plotFilesep.png
    def getExtNonCatalog(self, filterIn, midnight, starFile, plotFile=None,
				  amLbound=None, amUbound=None, magLbound=None, magUbound=None,
				  compareMags=False):
		#clear self.airmasses and self.magnitudes
		self.airmasses = []
		self.magnitudes = []
		
		#Compute night start time and night end time for night during which 'time' falls
		nightStartTime = midnight - dt.timedelta(hours=12)
		nightEndTime = midnight + dt.timedelta(hours=12)
            
		#Make a list, 'stars', of all stars at ras and decs listed in starFile
		f = open(starFile)
		line = f.readline()
		ra = []
		dec = []
		while line:
			word = line.split(' ')
			ra.append(float(word[0]))
			dec.append(float(word[1]))
			line=f.readline()
		stars = []
		for n in range(len(ra)):
			stars.append(self.database.getStar([ra[n], dec[n]]))
		
		#Make a list, 'usableStars', of all stars for which we have observations
		#on the input night and filter
		usableStars = []
		for star in stars:
			if (len(self.database.getStarExposures(star=star, filter=filterIn,
				startObsDate=nightStartTime, endObsDate=nightEndTime,
				minAM=amLbound, maxAM=amUbound, minMag=magLbound,
				maxMag=magUbound)) <> 0):
				usableStars.append(star)
		
		#Turn self.airmasses and self.magnitudes into lists of lists of size 
		#[number of stars][number of observations of star]
		enoughObservations = False
		for i in range(len(usableStars)):
			length = len(self.database.getStarExposures(
		                star=usableStars[i], filter=filterIn,
		                startObsDate=nightStartTime, endObsDate=nightEndTime,
		                minAM=amLbound, maxAM=amUbound, minMag=magLbound,
		                maxMag=magUbound))
			if length > 1: enoughObservations = True
			self.airmasses.append(numpy.zeros(length))
			self.magnitudes.append(numpy.zeros(length))
		if not enoughObservations: return ['Only one observation per star.  '+
			'Must have at least one star with multiple observations.']
		
		#Assign airmass and magnitude values to lists:
		#    airmasses[i][j] = airmass of jth exposure of ith star
		#    magnitudes[i][j] = instrumental magnitude of jth exposure of ith star
		allExpStars = list()
		i = 0
		j = 0
		while i < len(self.airmasses):
			j = 0
			expStars = list()
			expStars = self.database.getStarExposures(star = usableStars[i],
				filter=filterIn, startObsDate=nightStartTime,
                endObsDate=nightEndTime, minAM=amLbound, maxAM=amUbound,
                minMag=magLbound, maxMag=magUbound)
			while j < len(self.airmasses[i]):
				self.airmasses[i][j] = expStars[j].airmass
				self.magnitudes[i][j] = expStars[j].magnitude
				j += 1
			i += 1
		
		if compareMags: self.compareMagnitudes(allExpStars)
		
		#Compute various sums needed to calculate m0s and k
		mxIJ = 0
		mJ = 0
		xJ = 0
		x2J = 0
		x2IJ = 0
		N = []
		i = 0
		while i < len(self.airmasses):
			N.append(0)
			i += 1
		mx_N = 0
		x2_N = 0
		i = 0
		j = 0
		while i < len(self.airmasses):
			j = 0
			mJ = 0
			xJ = 0
			x2J = 0
			while j < len(self.airmasses[i]):
				mxIJ += self.magnitudes[i][j] * self.airmasses[i][j]
				x2IJ += self.airmasses[i][j]**2
				mJ += self.magnitudes[i][j]
				xJ += self.airmasses[i][j]
				x2J += self.airmasses[i][j]**2
				j += 1
			N[i] = len(self.airmasses[i])
			mx_N += (mJ*xJ)/float(N[i])
			x2_N += (xJ**2)/float(N[i])
			i += 1
		
		if x2IJ==0:
			return ['All observations made at airmass 0.  Cannot compute extinction.']
		print mxIJ, mx_N, x2IJ, x2_N
		numerator = mxIJ - mx_N
		denominator = x2IJ - x2_N
		self.k = numerator/denominator
		
		self.m0 = []
		mSum = 0
		xSum = 0
		i = 0
		while i < len(self.airmasses):
			j = 0
			mSum = 0
			xSum = 0
			while j < len(self.airmasses[i]):
				mSum += self.magnitudes[i][j]
				xSum += self.airmasses[i][j]
				j+=1
			self.m0.append((mSum - self.k*xSum)/N[i])
			i+=1
		
		#Compute correlation coefficient
		ssxx = 0
		ssyy = 0
		ssxy = 0
		xmean = 0
		ymean = 0
		for n in range(len(self.airmasses)):
			for m in range(len(self.airmasses[n])):
				xmean += self.airmasses[n][m]
				ymean += self.magnitudes[n][m]
				denominator += 1
		xmean /= denominator
		ymean /= denominator
		for n in range(len(self.airmasses)):
			for m in range(len(self.airmasses[n])):
				ssxx += (self.airmasses[n][m] - xmean)**2
				ssyy += (self.magnitudes[n][m] - ymean)**2
				ssxy += (self.airmasses[n][m] - xmean)*(self.magnitudes[n][m] - ymean)
		r2 = (ssxy**2)/(ssxx*ssyy)
		
		#Make plots
		if isinstance(plotFile,str):
			#Plot all observations and lines on one plot
			m = 0
			while m < len(self.airmasses):
				plt.plot(self.airmasses[m], self.magnitudes[m], 'ro')
				plt.plot([min(self.airmasses[m]), max(self.airmasses[m])],
					   [self.m0[m] + self.k*min(self.airmasses[m]), self.m0[m] +
					    self.k*max(self.airmasses[m])], 'b-')
				m+=1
			plt.suptitle('Non-catalog: k = '+str(self.k))
			plt.xlabel('airmass')
			plt.ylabel('m')
			plt.savefig(plotFile+'all.png')
			
			#Plot each star in a separate subplot
			minx = min([min(l) for l in self.airmasses])
			maxx = max([max(l) for l in self.airmasses])
			xrange = maxx-minx
			miny = min([min(l) for l in self.magnitudes])
			maxy = max([max(l) for l in self.magnitudes])
			yrange = maxy-miny
			fig = plt.figure()
			fig.set_figheight(len(self.airmasses))
			fig.subplots_adjust(hspace=1.5)
			sp = list()
			for n in range(len(usableStars)):
				am = []
				mag = []
				for m in range(len(self.airmasses[n])):
					am.append(self.airmasses[n][m])
					mag.append(self.magnitudes[n][m])
				sp.append(fig.add_subplot(len(usableStars),1,n+1,
					title=usableStars[n].objname+', m0 = '+str(self.m0[n]), ylabel='m_inst'))
				sp[-1].plot(am, mag, 'bo')
				sp[-1].plot([minx, maxx],
					[self.m0[n]+self.k*minx, self.m0[n]+self.k*maxx], 'r-')
				sp[-1].set_xlim(minx - xrange*.1, maxx + xrange*.1)
				sp[-1].set_ylim(miny - yrange*.1, maxy + yrange*.1)
			sp[-1].set_xlabel('airmass')
			plt.savefig(plotFile+'sep.png')
		
		#Create array consisting of m0.append(k)
		for n in range(len(self.m0)):
			toReturn={'m0'+str(n): m0[n]}
		toReturn['k': self.k]
		toReturn['corr. co.': r2]
		
		return toReturn
    
    def compareMagnitudes(self, expStars=None):
		f = open('compmag.txt', 'w')
		#Assume all exposure stars in same band
		band = expStars[0].exposure.filter
		i1 = self['filters'].index(band)
		i2 = self['filters'].index(self['referenceBand'])
		if i1 < i2:
			c1 = band
			c2 = self['referenceBand']
		elif i1 > i2:
			c1 = self['referenceBand']
			c2 = band
		else:
			c1 = self['filters'][i2-2]
			c2 = self['referenceBand']
		#Create a list of exposures
		exposures = list()
		for expStar in expStars:
			exp = expStar.exposure.filename
			hasExp = False
			for exposure in exposures:
				if exposure == exp: hasExp = True
			if not hasExp:
				exposures.append(exp)
		
		#Write comparisons:
		#Star1	Star2		m1-m01	m2-m02	(m1-m2)-(m01-m02)
		for exposure in exposures:
			f.write('Exposure: '+exposure+'\n')
			f.write('Star1\t\tStar2\t\tm1-m01\t\tm2-m02\t\t(m1-m2)-(m01-m02)\tcolor\n')
			for expStar1 in expStars:
				if expStar1.exposure.filename==exposure:
					for expStar2 in expStars:
						if (expStar2.exposure.filename==exposure and
							expStars.index(expStar1)!=expStars.index(expStar2)):
							m1 = expStar1.magnitude
							m01 = expStar1.star.magnitudes[band]
							m2 = expStar2.magnitude
							m02 = expStar2.star.magnitudes[band]
							color = 'n/a'
							if self['referenceBand'] in expStar1.star.magnitudes:
								color = expStar1.star.magnitudes[c1]-expStar1.star.magnitudes[c2]
							f.write(expStar1.star.objname+'\t'+
								  expStar2.star.objname+'\t'+
								  str(m1-m01)+'\t'+str(m2-m02)+'\t'+
								  str(m1-m2-m01+m02)+'\t\t'+str(color)+'\n')
			f.write('\n')
		
if __name__ == "__main__":
    print "This module is not intended to be run as a script."
